global KaiA KaiB KaiC m kxy0 kxyA Khalf k1 k2

%constants:
KaiA=1.3;                        % original concentration of KaiA
KaiB=3.4;                        % original concentration of KaiB
KaiC=3.4;                        % original concentration of KaiC
m=1;                             % parameter for adjustment
kxy0=[0 0 0 0.21;
      0.12 0 0.31 0;
      0 0 0 0.11;
      0 0 0 0];                   % matrix(4,4)
kxyA=[0 0.332923 0 0.0798462;
      0.123 0 -0.319385 0;
      0 0.505692 0 -0.133077;
      0.479077 0 0.0532308 0];    % matrix(4,4)
Khalf=0.43;                       % K_1/2 in equation
k1=60;                            % k_1 in the reaction of "S + KaiB = S-KaiB"
k2=40;                            % k_-1 in the reaction of "S + KaiB = S-KaiB"

% 1,2,3,4 represent T,ST,S,U respectively

tspan=[0,200];
y0=[0.68;
    1.36;
    0.34;
    KaiC-0.68-1.36-0.34;
    0.00];                         % original concentration: T=0.68,D=1.36,S=0.34,U=1.02,S-KaiB=0 , according to the paper

[t,Y]=ode45('odefunction_Amax',tspan,y0);

% plot the concentration of the four substantces
figure; plot(t,Y(:,1)*100/KaiC,'g',t,Y(:,2)*100/KaiC,'b',t,(Y(:,3)+Y(:,5))*100/KaiC,'r',t,Y(:,4)*100/KaiC,'k','LineWidth',2);
legend('T','ST','S','U'); xlabel('time / h','FontSize',14); ylabel('%KaiC','FontSize',14); axis([0 100 0 100]);

% proof of the rapid equilibrium
KK=KaiB*ones(length(Y(:,5)),1);
K2=k2*ones(length(Y(:,3)),1);
R=(Y(:,5)./Y(:,3))./((k1*(KK-Y(:,5)))./(K2+k1.*Y(:,3)));   
figure;plot(t,R);  % R is around 1 means that (Y(:,5)./Y(:,3)) can be seen as equal to ((k1*(KK-Y(:,5)))./(K2+k1.*Y(:,3)))
legend('R'); xlabel('time / h','FontSize',14);ylabel('ratio','FontSize',14); axis([0 100 0 2]);

% the fluctuation of KaiA
A=max(0,KaiA-2*m*Y(:,5));  % as in equation
AA=[A t];                  % create a variation to store A(t)
savefile='storeA.mat';
save(savefile,'AA');       % save A(t) in a .mat file for usage in Stochastic_ST.m and Stochastic_U.m
figure;plot(t,A,'LineWidth',2);
legend('A'); xlabel('time / h','FontSize',14);ylabel('concentration / \muM','FontSize',14); axis([0 100 0 0.7]);

Flow;   % plot the flow of the circle
